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Abstract 

We reformulate the Lanczos algorithm for quantum wave function propagation in terms of varia- 
tional principle. By including some basis states of previous time steps into the variational subspace, 
the resultant accuracy increases by several orders. Numerical errors of the alternative method ac- 
cumulate much slower than that of the original Lanczos method. There is almost no extra numeric 
cost for the gaining of the accuracy, i.e., the accuracy increase needs no extra operations of the 
Hamiltonian acting on state vectors, which are the major numeric cost for wave function propa- 
gation. A wave packet moving in a 2-dimensional Henon-Heiles model serves as an illustration. 
This method is suitable for small time step propagation of quantum wave functions in large scale 
time dependent calculations where the operations of the Hamiltonian acting on state vectors are 
expensive. 
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I. INTRODUCTION 



Propagation of quantum wave functions, i.e., direct integration of the time dependent 
Schrodinger equation, is a fundamental numeric task. This time dependent method exhibits 
many numeric advantages in first principle calculations. For example, one is able to extract 
energy spectrum efficiently from the correlation function via Fourier transformation Q], or 
more advanced filter diagonalization algorithm |2j. The efficiency is more evident when one 
needs excited energy spectrum in large scale first principle calculations. 

For wave function propagation, the most expansive numeric operations are products of 
the Hamiltonian matrix and state vectors, namely, the Hamiltonian operator acting on 
state vectors. It is a long standing efforts to develop efficient algorithm for wave function 
propagation that uses minimum number of such matrix-vector product operations. 

it operator method 



7, 81 is a robust and flexible scheme for wave 



Among the popular algorithms, such as split operator method 
expansion method lal, the Lanczos method 

H^—fl . . 

including the time dependent Hamiltonian [10( , and there is virtually no need for preparing 
knowledge about the considered system to apply the Lanczos method. Furthermore, the 
performance of the Lanczos method is relatively insensitive to the considered system and 
the initial wave function. In many cases, the Lanczos method is the best choice to do time 
dependent calculations. For example, if one need to compute a quantity (such as the entropy 
of a subsystem) changing continuously with time, and the corresponding Hamiltonian is not 
suitable to break into two parts to apply the split operator algorithm, one may consider to 
employ the Lanczos method for the task. 

The Lanczos algorithm transforms a Hermitian matrix into a tri-diagonal form itera- 



tively 



111 ]. It has many applications in first principle calculations, see, e.g. 
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The basic idea of wave function propagation by Lanczos method is to solve the Schrodinger 
equation for a given small time step in a low dimensional subspace, namely, the Krylov 
subspace. The basis states of the Krylov subspace, {^o, • ' ' ;VVt}, are generated by the 
Lanczos iteration, Hij)j = /^.i^j-i + + where = is the expectation 

value of the Hamiltonian H with respect to the vector /?$ is the norm of the vector 
Hipi — — aiipi with /3 = and if} being the wave function obtained from previous 

step. 
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The dimension of the Krylov subspace is usually less than 10 in most cases. This dimen- 
sion depends on the time step and the accuracy requirement. Higher accuracy needs either 
small time step or large Krylov subspace. For a given dimension of the Krylov subspace, 
the error accumulates linearly with time in Lanczos method 9]. If one needs wave function 
in longer time scale, one must increase accuracy of each time step to keep the error of the 
final wave function within required range. This means the numeric operations are not simply 
linearly proportional to the time. There is an optimal choice for the dimension of the Krylov 
subspace and the time step to reach the accuracy requirement of the final state. However, 
practical situations, e.g. calculations of correlation function, often need other time steps. 

The dimension of the Krylov subspace, or the number of matrix- vector product operations 
in a time step, is a key factor to affect the numeric cost for the Lanczos propagation scheme. 
Like other algorithms, the operations of the Hamiltonian acting on the state vectors are 
the major numerical cost of Lanczos propagation scheme. Higher accuracy demands for 
more such Hamiltonian operations. In the view point of efficiency, one should keep the 
number of the Hamiltonian operations as small as possible for a given time step and accuracy 
requirement. 

In this paper, we present an alternative to the Lanczos method. It can improve the 
accuracy by several orders without extra matrix- vector product operations. To this end, we 
reformulate the Lanczos method in terms of variational principle with the Krylov subspace 
being the variational subspace. The basic idea of improvement is to enlarge the variational 
subspace by including basis states of previous time steps into the variational subspace. Since 
the required matrix elements are already calculated in previous time steps, such enlargement 
of the variational subspace has virtually no extra numeric cost. In fact, including basis 
states of previous steps into the variational subspace is an efficient method for iteratively 
diagonalizing large matrix 



II. VARIATIONAL APPROACH TO THE LANCZOS PROPAGATION SCHEME 

We first note that one is able to formulate the Lanczos propagation scheme from the 
variational principle. For short time At, one can approximate the evolution operator U(t) = 
exp(—iH At / h) by a polynomial of the Hamiltonian operator H. In other words, one can 
approximate the wave function at time t + At, ip(t + At) = U(At)i/)(t), by a vector in the 
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Krylov subspace spanned by {H l ip{t), i = 0,1, • ■ ■ ,n — 1}, 

^(t + At) = ^c i H i V(t). (1) 

i=0 

The expansion coefficients q are yet to be determined by the variational principle, 

-^(t + At)|^|- - JET|V(* + At)) = 0. (2) 
OCi at 

One solves the resultant equation of motion to obtain the expansion coefficients. Instead of 
{H l if)(t)}, the basis states in original Lanczos scheme are generated by the Lanczos iteration, 
and the resultant Hamiltonian matrix in the variational subspace is tri-diagonal. Of course, 
the subspace generated by Lanczos iteration is the same as that spanned by states {H l if)(t)}. 

Another way of arriving at the ansatz (P) is to relate H l ip{t) with the z-th order derivatives 
of state ip(t) with respects to the time, ih-j^i/jft) = H l vp(t). One can approximate the state 
at time t + At, ip{t + At), by a linear combination of the state ip(t) and its derivatives up 
to order n. Such observation may be useful for time dependent systems. 

A direct way to improve accuracy of the Lanczos propagation scheme is to enlarge the 
variational subspace, i.e., adding more states into the basis states of variational subspace. 
This usually means more numerical operations of matrix-vector product to obtain the basis 
states and the correspondent matrix elements of the Hamiltonian. However, there exist a 
way to enlarge the variational subspace with virtually no extra numeric cost. 

We achieve this by making using of matrix-vector product operations of previous steps. 
To this end, we first note that one can propagate a wave function backward from if)(t) to 
ip(t — At). In other words, a wave function at time t — At, ip(t — At), is approximately 
a linear combination of states {H l ifj(t), i = 0, • • • , m}. Here m is the number of matrix- 
vector product operations in a time steps. At time t, we already have states H k ip(t — At), 
k — 0,1, ■ ■ ■ ,m. Each of these states is approximately equivalent to a linear combination of 
states {H k+i ip(t), i = 0, 1, • • • , m}. If one includes some of these states into the variational 
subspace, one has effectively products of higher order power of the Hamiltonian acting on 
the state vector ip{t). Using the same arguments, one can also include H k ijj(t - 2 At), 
H k ip(t — 3 At), ■ ■ • , into the variational subspace. 

Implementation of the above scheme is straightforward. It involves similar procedures 
as that of original Lanczos propagation scheme. Each step of propagation is to solve the 
Schrodinger equation in the variational subspace. The variational subspace is spanned by the 
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basis states of the original Krylov subspace, H l il)(t), and some basis states of previous time 
steps, H l tp(t — A; At). The solution of the Schrodinger equation in the variational subspace 
determines the expansion coefficients of next step's wave function with respect to the basis 
states of the variational subspace. The practical calculation includes the following steps: 

(1) Choose the time step At and the dimension of the variational subspace, n, as well as 
the number of matrix-vector product operations, m, for each step. The basis states of the 
variational subspace are 0^ = H l ip(t — kAt) with i — 0, 1, • • • , m, k — 0, 1, • • • , K. Here 
m < n, and ip(t) is the current state at time t. For practical applications, it is enough to set 
K = 1 and n < 2m, i.e., one usually needs only some of last step's basis states to form the 
variational subspace. In our implementation of choosing previous time steps' basis states, 
the order is H m ip(t-At), H™' 1 ^^- At), This is because that H m i)(t-At) is closer to 
H l ip(t), (i > m), than other states, and has less overlap with the basis states of the original 
Krylov subspace, H' l ^(t), (i < m). 

(2) Calculate the matrix-vector products (f$ = H l ip(t), % — 1, • • • , m + 1. These states 
and some states obtained in previous time steps form the basis states, 0^ = H l il)(t — kAt), 
of the variational subspace. Here 0o' n+1 ' 1 is only for calculation of the Hamiltonian's matrix 
elements in the variational subspace. Calculation of the m + 1 matrix-vector products in 
this step consumes the major CPU time of the whole procedure. 

(3) Calculate the matrix elements of the Hamiltonian in the variational subspace, Ttikji = 

= ($^|0! j+1) )> and the overlap between the basis states, S ik ji = (fc^l^). 
Here is normalized form of . For indices k > 0, / > 0, the matrix elements of H and 
S are already calculated in previous time steps, one needs only to calculate the terms Tiiojo, 
and S i0 jo in this step. The trade off of reusing previous steps' matrix element is that the 
basis states is not orthogonal with each other. 

(4) Solve the Schrodinger equation in the variational subspace 

ihS^C = HC (3) 
at 

to obtain the expansion coefficients C = (c^\ ■ ■ ■ , c£\ ■ ■ ■ ) T of next time step's wave function 
with respect to the basis states The computation cost in this step is negligibly small 
in comparison with other step's operations. 
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(5) Perform linear combination of the basis states to form next step's wave function 



At the first time step, t — 0, there is no previous basis states. All the basis states are 
formed by states H 1, i/j , % = 0, 1, • • • , n — 1, with ip being the initial state. In next step, 
we remove the first m + 1 states ip , Hip , ■ ■ ■ , H m t/j from the basis states, and add another 
m + 1 states ip(At), Hip (At), ■ ■ ■ , H m ^(At) into the basis states. In following time steps, 
we update the basis states of the variational subspace in the same way, i.e., replacing m + 1 
oldest basis states with states {H l ip(t), % = 0, • • • , to}. 

In case n — m + 1, i.e., the variational subspace including no previous time step's basis 
states, the above procedure is essentially the same as the original Lanczos propagation 
scheme. In such case, numerical cost, storage requirement, and resulted accuracy are indeed 
the same as the original one. In the original Lanczos scheme, the Hamiltonian matrix in 
the variational subspace is tri-diagonal, and the overlap matrix is unit. Since the dimension 
of the variational subspace is usually small, such difference results in virtually no extra 
numerical cost. Similar to the original Lanczos scheme, the above procedure is more suitable 
for small time step propagation which needs only small variational subspace. 

The storage requirement is also similar to the original Lanczos scheme. One needs to store 
the basis states of the variational subspace, as well as information about the Hamiltonian. 
Other informations, such as matrix elements of the Hamiltonian in the variational subspace 
and the overlap matrix, need little memory. 

For an given time step and accuracy requirement, there is an optimal choice of the di- 
mension, n, of the variational subspace and the number, m, of the matrix-vector product 
operations in a time step. If m is inadequately small, i.e., one includes too many previous 
time steps' basis states into the variational subspace, the overlap matrix may become sin- 
gular. This means that there is a limit accuracy for a given time step At and the number of 
matrix-vector product operations in a time step. We use this property of the overlap matrix 
to determine the dimension n for a given m and At. 




(4) 
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III. NUMERICAL RESULTS 



We test the performance of the alternative Lanczos method via Henon-Heiles model. It 
is a particle of unit ,„a S8 m ov,„ g in the 2-d lm e„ 8 io„al Henon-Heiles potential Q, „(*, ,) = 
\{u 2 x x 2 + u 2 y 2 ) + Xy(x 2 + r/y 2 ), where u x = 1.3, uo y = 0.7, A = —0.1, rj = 0.1, and the 
Planck constant is set to 1. This system has a chaotic classical limit. It is widely used to 
study the quantum-classical correspondence. Similar to Ref. |9j, we estimate the accuracy 
of the current method by the overlap between the numeric result and the "exact" result. We 
obtain the "exact" result via Chebyshev expansion method [5]. The Chebyshev method is 
a global propagator that can reach an accuracy of machine's limit with a single time step. 

In Figure 1, we show the auto-correlation function (if)(t)\ip(0)), i.e., overlap between 
initial state ip(0) and the state at time t, ip(t). The initial state is a Gaussian wave packet 
whose center positions are (2.0,2.0), and center momentums vanish. We use a 64 x 64 
grid to represent the 2-dimensional wave function in spatial representation. The action of 
momentum operator on the state is performed via Fast Fourier Transformation (FFT) to 
transform the state into momentum representation. Thus the action of the Hamiltonian 
operator on a state vector needs two FFTs to transform the state back and forth between 
coordinate and momentum representations. Such matrix-vector product operation is the 
major numerical cost of the wave function propagation. The thin solid line in Fig. 1 is 
a well converged result for comparison. The dashed line is result of alternative Lanczos 
method, and the thick solid line is result of the original Lanczos method. Both the original 
and alternative Lanczos methods use 2 matrix-vector product operations in one time step 
to obtain the auto-correlation function in Fig. 1. The time step is At = 0.02. It is evident 
that 2 matrix-vector product operations are not enough to converge for the original Lanczos 
method. In fact, one needs at least 4 to 5 matrix- vector product operations in a time step 
to make the original Lanczos method converge. On the other hand, the dashed line from 
alternative Lanczos method is almost indistinguishable from the well converged result. Here 
we include 8 previous basis states, {Wiplt — kAt), % — 0, 1, 2, k = 1, 2, 3}, into the variational 
subspace (ip(t — 3 At) is not included). The total dimension of the variational subspace is 
11. 

By reducing the time step to At = 0.01, one can obtain the above converged auto- 
correlation with only a single matrix-vector product operation in a time step. For such 
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time step, we achieve similar accuracy as that in Fig. 1 by including 5 previous time steps' 
basis states into the variational subspace. In contrast, for such time step and accuracy, the 
original Lanczos method still needs about 3 to 4 matrix-vector products in one time step. 

Generally, one can increase the accuracy with virtually no extra numeric cost by including 
more previous basis states into the variational subspace. However, in calculation of Fig. 1, 
including more than 8 previous basis states makes the overlap matrix singular, i.e., the basis 
states are no longer independent from each other. In fact, for a given time step At and 
the number m of matrix-vector product in a time step, there is always a limit number of 
previous basis states that one can include into the variational subspace. In other words, the 
time step At and the number m determine limit of the accuracy. 

In our test calculations of Fig. 1, the overlap matrix becomes singular occasionally, 
when this happens, one can simply remove the non-independent basis vectors. This can be 
done by, e.g., Cholesky decomposition of the overlap matrix. In our implementations, we 
replace each non-independent vector by one more state H l ip{t) of the Krylov subspace. Such 
treatment preserves the accuracy at the expense of one extra matrix- vector product. 

Practical implementation of the alternative Lanczos method is indeed more stable and 
robust than the calculation of Fig. 1. In fact, calculation of Fig. 1 includes several previous 
time steps' basis states into the variational subspace. This treatment almost reaches the 
accuracy limit with 2 matrix-vector product operations in a time step. Practically, one 
needs only including some of last step's basis states into the variational subspace. This can 
increase the accuracy by several orders with virtually no extra numeric cost. The overlap 
matrix, and thus the whole numeric procedures, are usually well behaved. 

Fig. 2 shows the numeric errors accumulate with time for situations similar to practical 
calculations. The solid lines are results of the alternative Lanczos method, and the dashed 
lines are results of corresponding original Lanczos method which uses the same number of 
matrix- vector product operations. Same as Ref. [j]], we use the overlap between numeric 
result ipnumeric and "exact" result ipexact as the measure of error, 



We obtain the "exact" state vector by Chebyshev expansion method with accuracy of ma- 
chine's limit. We propagate the "exact" states with a time step AT = 4, and using 1024 
Chebyshev polynomials to expand the evolution operator exp(— iHAT/h). This is well 
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beyond the accuracy requirement of the machine's limit. From our tests, 512 Chebyshev 
polynomials are indeed well converged. From top to bottom of fig. 2, N is the total dimen- 
sion of the variational subspace for the alternative Lanczos method, and m is the number 
of matrix-vector product operations in a time step for both methods. The time step is set 
to At = 0.02 for both methods. The initial states, as well as states representation are the 
same as that of Fig. 1. It is easy to see that the behavior of the original Lanczos method is 
similar to that described in Ref. i.e., the error accumulates about linearly with the time 
t, and the accuracy increases quickly with the number m. 

It is evident that, when including some of last step's basis states into the variational 
subspace, the accuracy improves drastically. We see that, for m = 5, the alternative method 
is about 5 orders more accurate than the original one after time t « 10 3 . And for m = 6, 
the alternative method is about 4 orders more accurate after time t ~ 5 x 10 3 . The original 
method is well converged within the time scale t < 10 4 for m = 7. Even so, the alternative 
method is still about one order more accurate after time t ~ 2 x 10 4 . Another encouraging 
property of the alternative method is that its error accumulates much slower than the original 
one. This means numeric result of the alternative method is more reliable in long time scale. 

From Fig. 2 and other test calculations, we conclude that the alternative Lanczos method 
is suitable for small step wave function propagation. It improves the accuracy by several 
orders with almost no extra numeric cost. In calculation of Fig. 2, the variational subspace 
is about 10 dimensional, and includes only 4 last step's basis states. For such setting, the 
resultant overlap matrix and thus the overall numeric procedure are well behaved. In fact, 
the included basis states, H l ip(t — At), from previous step play the role of High order power 
of the Hamiltonian acting on the state vector, H l ip(t), of the original Lanczos method. From 
the dashed lines in Fig. 2, we see that basis states, H l Tp(t) (i — 0, 1, • • • , m), of the original 
Lanczos method span the major part (> 99%) of the exact wave function ip(t + At). Other 
basis states, H l ip{t — At), included from previous step span only very small portion (< 1%) 
of the exact wave function. Thus, these included basis states H %r i\}{t — At) have a relatively 
lower accuracy requirement to represent the high order terms H l ip(t), (i > m). This explains 
the success of the alternative method. 

In general, implementation of the alternative method is stable and robust, provided that 
the basis states {H l ip(t)} span major part of the next step's wave function if)(t + At). In fact, 
for a given time step At, the accuracy of both alternative and original Lanczos methods is 
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determined by the dimension of the variational subspace n. One chooses At and n in the same 
way as that of the original Lanczos algorithm. In the alternative implementation, one must 
specify additionally the number, m, of matrix-vector production in a time step. It is usually 
enough to set m larger than half of n, i.e., n/2 < m < n. If m is too small, and one includes 
too many basis states, {H l ip(t — kAt)}, of previous steps into the variational subspace, 
these basis states may be not independent. Even this happens, the alternative method still 
works. If a basis state from previous step is linearly dependent on other basis states of the 
variational subspace, we replace this state by an extra state H l ip(t) with i > m. This keeps 
the dimension of the variational subspace, and hence the accuracy of the implementation. 
The expense of such treatment is one more matrix- vector product operation for each linearly 
dependent state. We implement this treatment during the Cholesky decomposition of the 
matrix S which is a necessary step to solve Eq. (J3J). When the basis states are linearly 
dependent, the overlap matrix S becomes singular. The Cholesky decomposition of S can 
find all the non-independent states. By properly specify the threshold value of the Cholesky 
decomposition, this numerically cheap operation can even find the states that are close to 
linear superposition of other basis states. 



IV. CONCLUSIONS 



In summary, we present an alternative to the Lanczos method for quantum wave func- 
tion propagation in terms of variational principle. This method approximates short time 
evolution operator, U(t) = exp(—iHAt/h), by a polynomial of the Hamiltonian. In other 
words, the wave function i[)(t + At), resulted from a small time step propagation from 
ip{t), ip(t + At) = U(At)ip(t), is approximately a vector in the Krylov subspace spanned 
by {H l ip(t),i = 0, • • • — One can employ the variational principle to determine the 
expansion coefficients. The original Lanczos method needs to calculate all the basis states 
H l %p{t) explicitly. Construction of theses basis states is the major numeric cost. The al- 
ternative method needs only to calculate some of the basis states, H l ip{t), % = 0, • • • ,m, 
which span the major part of the wave function ip(t + At). We use basis states of previous 
step to play the role of the other basis states, H l i/j(t), i > m. Practically, it is enough to 
include some of last step's basis states, H l ip{t — At), i < m, into the variational subspace. 
The accuracy of the alternative method is several orders higher than the original Lanczos 
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method with same matrix-vector product operations in a time step. 

This alternative method is especially efficient for small time step wave function prop- 
agation. The error accumulation in the alternative method is much slower than that in 
the original Lanczos method, which increases about linearly with time. The efficiency of 
the alternative method comes from the fact that the basis states included from previous 
steps only span very small portion of the wave function, and thus the accuracy requirement 
for construction of these basis sates is relatively lower. This alternative method has useful 
applications in large scale time dependent calculations in which the numeric cost for the 
Hamiltonian acting on state vectors is expensive. 




20 40 t 60 80 

FIG. 1: Real part of auto-correlation function changes with the time. The thin solid line is a 
well converged result. The dashed and thick solid lines are results of the alternative and original 
Lanczos methods respectively. Each time step uses 2 matrix-vector product operations for both 
dashed and thick solid lines. 
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FIG. 2: Logarithms, base 10, of numeric errors, Err, accumulate with the logarithmic time t. 
Solid and dashed lines are the alternative and original Lanczos methods, respectively, m denotes 
the number of matrix-vector product operations in a time step for both methods, and N is the 
dimension of the variational subspace for the alternative method. 
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